arXiv:1505.03422vl [math.NA] 13 May 2015 


Mimetic spectral element method for Hamiltonian systems 


Artur Palha®’’*, Marc Gerritsma*’ 

^Eindhoven University of Technology, Department of Mechanical Engineering, section of Control Systems Technology P.O. 

Box 513 5600 MB Eindhoven, The Netherlands 

^Delft University of Technology, Faculty of Aerospace Engineering, Aerodynamics Group P.O. Box 5058, 2600 GB Delft, 

The Netherlands 


Abstract 

There is a growing interest in the conservation of invariants when numerically solving a system of ordinary 
differential equations. Methods that exactly preserve these quantities in time are known as geometric 
integrators. In this paper we apply the recently developed mimetic framework |461 152) to the solution 
of a system of first order ordinary differential equations. Depending on the discrete Hodge-* employed, 
two classes of arbitrary order time integrators are derived. It is shown that the one based on a canonical 
Hodge-* results in a symplectic integrator, whereas the one based on a Galerkin Hodge-* results in an energy 
preserving integrator. A set of numerical tests confirms these theoretical results. 

Keywords: Energy-preserving methods; Symplectic methods; Mimetic methods; Spectral element method; 
Collocation methods; Geometric integration 


1. INTRODUCTION 


In the last decades there has been a significant development of discretization methods for the solution of 
ordinary differential equations {ODEs). These equations are ubiquitous, appearing in particle accelerators, 
molecular dynamics, vortex particle methods, stock market modeling, celestial dynamics, etc. From a 
mathematical point of view, these problems can be generalized as a system of first order ODEs: 

^=h{y,t), y{to)=y°, y G C K” and t G [toP/] C K . (I.I) 


Where h{y, t) G M" and the region Ki is the phase space of the system. For the case of Hamiltonian systems, 

y = (p,g), q e 


q e 


and (1.1) takes the form: 


dt 


J ^\/H{y), y(to) = y G Ad C and tG/CM with J 



I 


m 


0 


( 1 . 2 ) 


where 1^ is the m x m identity matrix and h = J“^ViiI(y). 

Traditionally, the methods developed for the discrete solution of these problems derive from finite dif¬ 
ference principles and Taylor series expansions. This approach focuses on minimizing the local truncation 
error associated with the discretization. No matter how small this truncation error is, it is always finite, 
which might lead to a continuously growing global error. Over a sufficiently long time interval, the global 
error will inevitably reach unacceptable values and eventually the numerical solution will no longer resemble 
the analytical one. Due to this intrinsic local truncation error of the discretization process, one wishes to 
develop methods that are able to preserve the same qualitative properties as the continuous system. This 
will enable acceptable long time integration. 

By qualitative properties one refers, for example, to the undermentioned properties, |16j : 
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GeomeTBie staicittiiKDB: equations preserve geometric properties, such as orthogonality of lines or surfaces, 
volumes in phase space, etc. 

ConserWMdmd’aralsdng in time, certain well known quantities as mass, momentum, energy, or other less 
known quantities such as the radius of the trajectory or functions of the dynamical variables are 
preserved. 

SymmcMtan^^ systems remain invariant under specific transformations, e.g.: Galilean transformations 
(translations, reflexions and rotations), reversal symmetries (of which time reversal is an example), 
scaling symmetries (invariance under rescalings in either time or space). 

The key word to retain from qualitative properties is conservation. When we refer to qualitative properties 
we mean the invariance of certain quantities or geometric relations, during time evolution. These qualitative 
properties are explicit constraints that restrict the universe of allowable solution trajectories, discarding 
solutions which are not physically compatible with the system of equations one wishes to solve. 

A simple yet illustrative example for the relevance of geometric integration is presented, among others, 
in m for the discretisation of the harmonic oscillator with a forward Euler scheme. 

Example 1. Consider the two coupled ODEs that describe the evolution of an harmonic oscillator: 



dp 

dt 


= -q 


with q[to) = go, p(to) = Po and t S [toCf] ■ 


(1.3) 


It is possible to verify that all solutions are periodic and bounded and that [q^ +P^) is a conserved quantity of 
the evolution. If a forward Euler discretization is applied to the system, the following discrete time evolution 
expressions are obtained: 


{ qn+l = qn + ^tpn 
Pn+l =Pn- I^tqn 


with g(<o) = Qo, p{to) = Po ■ 


(1.4) 


It is straightforward to see that with this time stepping scheme the quantity [q^ + p^) grows geometrically 
at each time step: 

Qn+l +Pn+1 = (1 + + Pn) ■ 

Due to this feature of the forward Euler discretization, the discrete solutions of the harmonic oscillator will 
no longer be periodic and bounded. Moreover, since up to a multiplicative constant i^q^ + p^) is the energy 
of the system, conservation of energy is no longer satisfied. Another relevant point is that this discretization 
renders the discrete system time-irreversible, in contrast to the time-reversibility of the continuous system. 
This can be seen by using (I. 4 ) with a negative time step starting from (x„+i, g„+i)/ 


{ qn = g„+i - Atp„+i = + Atpn - At{pn - Atqn) = gn(l + At^) 

Pn = Pn+l + Atqn+1 = Pn - At{qn + Atp„) + Atqn+1 = Pn(l + At^) . 


(1.5) 


Where we have used { I. 4 ) to replace qn+i and pn+i. 
reversible. 

On the other hand, if the implicit midpoint rule is 
evolution expressions become: 


Since 7 ^ qn and pn ^ Pn, this method is not 
applied to this system, the alternative discrete 


time 

time 


Qn+l — gn 


4- At^ 
4 + At2 
4- At^ 
4 + At2 


■ Vn 


AAt 
4 +At2 
AAt 


Qn ■ 


Pn+l — Pn 


4 +At2 
2 


with q{to) = qo, p(to) = Po ■ 


( 1 . 6 ) 








For this discretization, the quantity is conserved, regardless of the time step used: 

9n+l Pn+1 Qn Pn ■ 

Contrary to the forward Euler discretization, the midpoint rule is capable of preserving the same qualitative 
properties of the continuous system, namely: the discrete solutions are periodic and bounded and the quantity 
(<?^ + P^) 'is conserved. Additionally, it is possible to show that the discrete system is time reversible, just 
like the original dynamics. 

This simple example expresses the relevance of judiciously selecting a numerical integrator. If correctly 
chosen, it is capable of preserving the qualitative properties of the continuous system. The approach that 
incorporates global characteristics of the problem being solved to construct discretization schemes with the 
same properties of the continuous one, gives rise to a class of time integrators referred in literature as 
geometric time integrators. 

In this paper, the focus is on autonomous Hamiltonian problems, i.e. ones such that h{y, t) = H (y). 
The scalar function H (y) is the Hamiltonian of the problem and its value is independent of time, 

H [y{t)) = H {yo), Vf > to • (1-7) 

In the case of isolated mechanical systems, the Hamiltonian and energy are associated. Another characteristic 
of Hamiltonian systems is that they are symplectic, see m for example. The exact flow map, ipr, associated 
to a system of ODEs is defined as (firiyo) ■= v{t), where yir) is the solution of the initial value problem at 
time t = T. A system of ODEs is symplectic if the associated Jacobian of the flow map, satisfies: 

Fripf Jp'riy) = J- ( 1 - 8 ) 

Symplecticity implies, for instance, that the flow map preserves phase space volume, which imposes strong 
constraints on the admissible solutions. This can be stated mathematically as: 

dg(<) A dp(t) = d9(to) A dp(to) , (1-9) 


where t > tg. 

If we apply equation (1.9) to the Euler discretization, (1.4), used in Example[^we see that: 
dg„+i A dpn+i = {dqn + Atdpn) A (dp„ - Atdg„) = dqn A dp„ + (dg„ A dp„) , 


showing that the Euler discretization is not symplectic. On the other hand, if we apply equation (1.9) to 
the midpoint rule discretization, (1.6), used in the same example we obtain: 

4At2 


4 — At^ 

dq„+i A dp„+i = -r—rTzdq„ A dp„ 
4 + 


4 +At2 


dqn A dpn = dqn A dp„ 


confirming the symplecticity of the midpoint rule. 

Since Hamiltonian systems are fundamental in the study of many physical systems, several discrete time 
integrators have been developed. As said before, energy conservation and symplecticity are key properties 
of Hamiltonian dynamics to satisfy at a discrete level. Ideally, a numerical integrator with discrete flow 
map (fAt such that y„^i = (pAt iVn) should be symplectic and exact energy preserving. Although highly 
desirable, it is very difficult to have a numerical integrator that preserves all geometrical properties. In this 
work we take as starting point the mimetic framework |46[ 152) and show that two numerical integrators of 
arbitrary order, one symplectic and the other exactly energy preserving, can be derived. 


1.1. Literature review 

Structure-preserving methods that aim to preserve at a discrete level the geometrical, topological and 
homological properties of the systems of equations being solved have gained increased popularity in the last 
decades, see m for a brief overview. Two main areas of development have emerged: one focussing on spatial 
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discretization (mimetic/compatible discretization) and another focussing on time discretization (geometric 
integration). 

Regarding spatial discretisations, mimetic/compatible discretizations have been proposed by several 
authors for the numerical solution of partial differential equations. Tonti’s work |69] was one of the first 
to establish the relation between differential geometry and algebraic topology in physical theories and its 
implications on the development of structure preserving discretizations. Tonti employed differential forms 
and cochains as the building blocks of his method. The relation between differential forms and cochains was 
established by the Whitney map (fc-cochains —>■ fc-forms) and the de Rham map (fc-forms —?> fc-cochains). 
These operators were initially introduced by Whitney in |72] and de Rham in [22) . respectively. The 
linear interpolation of cochains to differential forms on a triangular grid was established in |72j employing 
what is now known as Whitney forms. On quadrilaterals, Robidoux, |62j . and Gerritsma, EH, derived 
the edge basis functions, capable of arbitrary order interpolation and histopolation of differential forms. 
Subsequently, among others, Robidoux, [60l[6T], Hyman, [40ti42] . Steinberg, [66l[67], Shashkov, [64], Brezzi, 
[laila], Desbrun, [MlEaiSI], Lipnikov, [48], Perot, |55|, and Pavlov, [53], have proposed approaches in a 
finite difference/volume context. We highlight the ‘Japanese papers’ by Bossavit, (TUI HI], which serve as 
an excellent introduction and motivation for the use of differential forms in the description of physics and 
its use in numerical modelling. In a series of papers by Arnold, Falk and Winther, EHl], a finite element 
exterior calculus framework is developed. Bochev, |8], presents general principles of mimetic discretizations, 
common to finite element, finite volume and finite difference discretizations. Finally, Kreeft et ah, [46], and 
Palha et ah, [^, introduced a framework merging the ideas of Tonti and Bossavit, resembling a hybrid 
finite element/finite volume formulation of arbitrary order on curvilinear elements. 

On the topic of time discretization, namely the numerical integration of Hamiltonian systems, two main 
paths have emerged: symplectic and exact energy-conserving integrators. For a detailed overview of these 
methods the reader is referred to m- Symplectic integrators are widely known and can be traced back to 
the work of Vogelaere [23] and have been extended to Runge-Kutta methods by several authors, for example 
[231 IS31 EH]- Exact energy-preserving methods can be found in the works of Gonzalez (for Hamiltonian 
systems with symmetry), [33], Quispel et ah, [S3], and McLachlan et ah, [SD], (for general Hamiltonians). 
This class of methods, named discrete gradient methods, defines a discrete counterpart of the gradient 
operator ensuring that energy conservation is guaranteed in each time step. Projection methods are an 
alternative to discrete gradient methods. Hairer makes an extensive introduction to projection methods 
in [38] and later extends this method to include time reversibility in |35j . Unlike discrete gradient and 
projection methods, the averaged vector field method (AVF method) introduced by Quispel, [57], does not 
require the knowledge of the functional form of the conserved integral, requiring only information on the 
vector field. An extension of AVF methods to Runge-Kutta methods is presented in [18]. Subsequently, a 
modification of collocation methods has been introduced in |36j . extending the AVF method to arbitrary 
order. Following a different route, the time finite element method, initially introduced in [21[3nil3S] , was 
solved by a mixed formulation in [5] and extended to satisfy energy conservation in [7]. More recently, line 
integral methods, |14j , (of which Hamiltonian boundary value methods, |14L 115] , are a subclass) have been 
introduced. These methods exactly preserve the energy of arbitrary order polynomial Hamiltonian systems. 

With respect to discretization methods that combine both space and time it is relevant to mention Box 
schemes or Keller box schemes, which are a class of numerical methods originally introduced by Wendroff, 
for hyperbolic problems and later popularized by Keller, [331111], for parabolic problems. These methods 
are face-based and space and time are coupled by introducing a space-time control volume. They are known 
to be physically accurate and several successful applications have been reported, [201 mi m EH]- More 
recently, this method has been shown to be multisymplectic (Ascher, [6], and Frank, [28]). Also, Perot, [54] . 
established a relation between box schemes and discrete calculus and generalized it to arbitrary meshes. 

1.2. Outline 

Despite the increasing popularity and success of mimetic/compatible methods their extension to the 
solution of ODEs has received little attention and mainly concerns time-marching schemes for PDEs, [min]. 
Therefore, in this work we present how the mimetic framework, [461 152] . can be applied to the numerical 
solution of systems of ODEs. As seen in [331 [S2], all approximations lie in the discretisation of the Hodge-* 
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operator, which may be defined in different ways. In this work we show that a discrete Hodge-* obtained 
from a canonical Hodge-* operator results in a symplectic time integrator, whereas a Galerkin Hodge-* 
operator gives rise to an exact energy preserving integrator. 

The outline of the paper is as follows. In Section a brief introduction to differential geometry is 
presented. The focus will be on vector fields, differential forms and the operators that act upon them. In 
Section we introduce the required elements of algebraic topology, showing a clear parallel structure to 
differential geometry. In Section we establish a connection between the continuous world and the discrete 
by introducing the mimetic projection operator and the reconstruction of fc-forms. In Sectionthe mimetic 
integrators and their properties are discussed and in Section]^ the time integrators are applied to test cases. 
Finally, in Section the summary of this work and further applications are discussed. 


2. A brief introduction to differential geometry 

As mentioned previously, the objective of this work is to apply the mimetic framework developed in 
[461152| to the solution of systems of ordinary differential equations. The mimetic framework is based on 
differential geometry for the representation of the physical quantities and the operators that act on them. 
For this reason, a brief introduction to the key concepts needed to understand this approach will be given. 
For a more detailed presentation of this topic see, for example, [3 [HI HZl US]. Since time, t, is the only 
independent variable in the problems dealt with in this work, we give special attention to manifolds of 
dimension n = I. For this reason, each new concept is first introduced for manifolds of arbitrary dimension 
and then the particular case of n = I is presented. Throughout the text we will consider Riemannian 
manifolds. A/", of arbitrary dimension n and temporal Riemannian manifolds, T, of dimension one. 

Consider a smooth curve 7 ( 3 ) in M parametrized by s G [—cr, cr], cr > 0, such that 7 ( 0 ) = P £ Af. The 
derivative 7 ( 0 ) is a tangent vector at the point P € Af. For n-dimensional manifolds it is possible to find a 
set of n linearly independent tangent vectors at the point P, { ej |p ,..., e„|p}, see for example PQHT). This 
set spans a linear vector space called tangent space at the point P and denoted by TpAf. Any vector at P 
can be written as a unique linear combination of these basis vectors, i.e., 

n 

^Ip = ■ (2.1) 

i 


Where the coefficients u* are associated to the particular basis elements ei|p. If in the vicinity of the point 
P we define a local coordinate system, {x^,... ,a;”}, we specify the basis of the tangent space and in this 
case the basis vectors are generally denoted by ^ | 

( 2 . 1 ) for a local coordinate system as: 


and called coordinate basis vectors. We can then write 




_d_ 

dx^ 


It is possible to smoothly define a tangent vector at each point P G Af. This construction generates vector 
fields. To simplify the notation, in what follows we will suppress the explicit reference to the point P. 


For a temporal one dimensional Riemannian manifold, T, the tangent space at each time instant r, 
Tt-T, has dimension one and can therefore be spanned by a single coordinate basis vector {^}. Any 
tangent vector, v G T^T, can be represented as: 

d 

^ = ^df 

where the coefficient v is associated to the coordinate basis. 
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For a Riemannian manifold Af, there exists a metric gp, or simply g, that assigns to every point P € Af 
an inner product g = (•, ■)p on TpAf which depends smoothly on P, i.e. 

g : TpJ\f X TpJ\f —>■ M . 


For the coordinate basis vectors we define the metric coefficients as: 

Due to the bilinearity, the inner product between any two real tangent vectors v = J2i ^ s-nd u = 
can be written as: 


{v,u) 





■E 




dx^ ’ dxt 




(2.3) 

(2.4) 

(2.5) 


If we consider a one dimensional Riemannian manifold, T, there exists only one metric coefficient: 


511 


d d 
dt’ dt 


( 2 . 6 ) 


To simplify the notation, instead of gn we will use g to denote the metric coefficient for a one dimensional 
manifold. 


With any linear vector space V we can associate the dual space V* of linear functionals acting on V, 

Va€V*, a:V^R. 

In the same way, with TpAf we can associate the dual space TpAf of linear functionals acting on the tangent 
space: 

: TpAf M. 

We say that G TpAf, where TpAf is the cotangent space. The elements of the cotangent space are 
called covectors or alternatively exterior 1-forms. As was done for the tangent space, for a particular local 
coordinate system we can define a basis {dx^,... ,dx^} and represent these exterior 1-forms as: 

= aidx^ -I- • • • -I- a„dx" . 


These basis 1-forms are such that: 



where d* is the Kronecker-d. 

A differential 1-form, a^, is a smooth assignment of an exterior 1-form, 

ot^ ^ (xi,..., xffi Oi\(^Xi ,..., Xji^dx ‘cxji(^X\,..., Xjfjdx G TpAf , 

to each point, P G Af. We write G Af{Af) or simply G A^. 
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For a one dimensional manifold r, with TtT we can associate the dual space of linear functionals 
acting on the tangent space. In this case the basis consists of only one element {dt} and the exterior 
1 -forms are represented as: 

= adt , (2.7) 

where a is the t-component of the exterior 1-form. In this coordinate system we have that: 


dt 



= 1 . 


And therefore: 


(w) 



= av. 


The exterior product or wedge product, A, allows the construction of higher rank A:-forms, 1 < k < n, 
from 1-forms. Let Af he a manifold of dimension n, A* (Af) and A* (AA) be the space of fc-forms and /-forms, 
respectively, with k + I < n, then the wedge product. A, is a mapping: 

A : A'^ (Af) X A' (Af) A'=+' (Af), k + l<n 
which satisfies the following properties: 


A ^(m) = Q,(fc) A j(m) pd) ^(m) (Distributivity) (2.8a) 

A A 7 ^"*^ = A A A A 7 ^™^ (Associativity) (2.8b) 

A A = a A (Multiplication by functions) (2.8c) 

Q,(fc) /\ ^(0 = (— 1 )^* 13 d) /\ Q,(fe) (Skew symmetry) (2.8d) 

where € A^^^ (A/”), € A^*^ (AA) and € A^™) (A/”). A differential /c-form, € A^, is then an 

alternating /c-tensor: 

: TpAA X • • • X TpA/" ^ M , 

k times 


and we can write: 


Afc) _ 


= ^ ... ,a;„) da:*i A • • • A da;**" , 1 < ti < • • • < ifc < n . 


(2.9) 


The space of differential 0-forms, A(°^(A/'), is simply the space of smooth functions and we write: 

/3(°)(xi, ...,Xn)= Pixi, ... ,x„). (2.10) 


For the particular case n = 1 the wedge product takes the simpler form: 

a(i)A/3(°) =/3(0)Aa« =a/3dt, /^(o) Ay^^) = 7 ^°) A/3(°) = ( 7 / 3 )(°) and Aa^ = Aa« = 0 , 

where € A°(T) and € A^(T). 

Differential fc-forms naturally integrate over /c-dimensional manifolds. Let € A* (AA) and A4 C AA C 
M", with k = dim (A4), n = dim {AA) and k < n, 
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( 2 . 11 ) 


(aW,A4) := f 

JATk 

this represents a metric-free duality pairing. 


For the case n = 1 considered here we have: 

{/3^°\To) ■■= [ and (a^^^Ti) := [ . 

JTo JTi 

Note that To is a 0-dimensional manifold, therefore a point. By definition, the integral of a scalar 
function over a point is just the evaluation of the function at the point: = PiTo)- 


For Riemannian manifolds, Af, we can define a point-wise inner product of fc-forms, (•, •): 

(•,•) : A'=(Af) X A° (Af) . 


If and are two fc-forms expressed in the form (2.9), then the point-wise inner product is given by 




Where := (da:*,da:f) are the components of the inverse of the metric tensor, (2.2). 


For 1-dimensional manifolds, T, the point-wise inner products between 0-forms, G A°(T), 

written in the form (2.10) is given by 


=aP 


( 2 . 12 ) 


and the inner product between 1-forms, € A^(T), written in the form (2.7) is given by 

= aP - . 


(2.13) 


Note that we used g := gu, with gu as in (2.6). Therefore g^^ “ ^ ~ s' 


The Hodge-* operator is a mapping: 

* : A'=(Af) ^ A”-i(Af), 

that is induced by a combination of the wedge product and the inner product, 

A */?('=) := cr(") . 


(2.14) 


Where cr^"^ is the n-form, called volume form, such that cr^"^ is the volume of the manifold AT and 
*1 = 

The inner product, (•, •)i2(_A/), is defined as a mapping, 

(•,.)z.2(aa):A'=(AT)xA'=(AT)^M 


such that: 



L2(N) 


JN 



A ■ 


(2.15) 
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The Hodge-* operator, for the case of one dimensional Riemannian manifolds, T, takes the form: 

1 


* 1 = dt and * df = 


V9 


The inner product between 0-forms, € A°(T), written in the form (2.10) becomes 


and the inner product between 1 -forms, G A^(7~), written in the form (2.7) is 




Jj- 


dt. 


(2.16) 


(2.17) 


(2.18) 


Differentiation in differential geometry is represented by the exterior derivative, d, which is defined as a 
mapping, 

d: 

that satisfies the generalized Stokes Theorem: 



(2.19) 


Where AC is a (fc -I- l)-dimensional manifold and dN is its fc-dimensional boundary manifold, and G 
hA{dN). This equation contains, as special cases, the Newton-Leibniz, Stokes and Gauss theorems and hence 
generalizes the vector calculus V, Vx and V- operators to arbitrary manifolds. Moreover, this operator 
is completely topological (does not depend on the metric) and satisfies the Leibniz rule d(Q;(^^ A = 
da^^^ A -I- (— l)^'a^^^ A d/3^^^. An important point to note is that this theorem can be expressed using the 


duality pairing, ( 2 . 11 ), as: 


(da('=),AC) = {a''^\dN) . 


( 2 . 20 ) 


This duality property will play a fundamental role in the construction of the discrete exterior derivative. 


For 1-dimensional manifolds, the exterior derivative, d, is a mapping d : A*^ i—> A^, such that on a 
coordinate system it is defined as: 

da(o) := — dt, e . ( 2 . 21 ) 

at 

In this case, since n = 1, we have that da*^^^ = 0, G A^(T). Higher order differentiation can be 

defined by combining the exterior derivative and the Hodge-* operators, for example: 


d * da(°) 


da , , 1 da d 1 da , 

dt = d— — = — — — dt 
dt dt dt dt 


In this way one can construct a double de Rham complex associated to differential forms in 
1-dimensional Riemannian manifolds, (2.22): 



* Ai 

/K 

-* 

-AO 


( 2 . 22 ) 


9 

























There are two fundamental results to extract from the differential geometric formulation. First, from the 
generalized Stokes’ theorem in one dimension, it is possible to recognize that 1-forms are intrinsically related 
to integrals over time intervals whereas 0-forms are associated to evaluations at time instants. Moreover 
exact differentiation is possible if the points where 0 -forms are evaluated are the boundary of the line 
segments where 1-forms are integrated. Second, from the double de Rham complex, (2.22), one can see that 
there exist two distinct complexes, the top one and the bottom one, both being related to each other by the 
Hodge-* operator. It states, for example, that is associated to time intervals and that * 0 ^^^ is associated 
to time instants. Conventionally, these two quantities are considered the same quantity and are discretized 
in the same way. It is the authors’ opinion that this formal distinction is essential in the construction of 
a numerical discretization since, at the discrete level, the Hodge-* cannot be performed exactly. Therefore 
one needs to make use of dual grids to represent this dual complex. In this way all approximations lie in 
the discretisation of the Hodge-*, see [Ml SSI HU [521 Eg ES]. 


3 . A brief introduction to algebraic topology 

Following the mimetic framework presented in |4bL I52j , consider a temporal 1-dimensional domain, T, 
and its primal, D, and dual, D, grids, see Figure]^ 

Ri).! Ril.a Ri),(Af-i) Hi),n 


D 




Ro),o Ro),i 

Ro).i Ro).2 


T(0),{N-1) ''■(0),Af 

Ro),(Af-i) Ro),Ar 


D 


Ri).o Ri).i 


Rl).(iV-l) T(l),Ar 


Figure 1: Example of a primal, D, and dual, D, grids covering a 1-dimensional manifold. Notice that the number of nodes, 
I'he primal mesh is equal to the number of edges (line segments), in the dual mesh, and vice-versa. 


These grids consist not only of the points, 

{ro),z :* = 0, ...,V} and {f(o),* : i = 1, ■ ■., A} , 

as is common in many numerical discretizations, but also of the edges (line segments) connecting them, 

{ri),^ : * = 1, ■ ■ •, A} and : i = 0,..., A} . 

One important point to note is that the number of nodes, Rq) j, in the primal mesh is identical to the number 
of edges, ■?(!),i, in the dual mesh, and vice-versa. 

The fc-dimensional objects in D and D are called k-cells and we represent them by T(^k),i where k denotes 
the dimension of the object (in our case fc = 0,1) and i is a label that distinguishes different objects. 

It is possible to inner orient a time interval (a line) by defining a direction that goes from the preceding 
time instant to the following one (as in the natural time sequence), see top left of the Inner orientation block 
in Figure]^ One can also inner orient a time interval in the opposite direction (as in time reversal), see top 
right of the Inner orientation block in Figure]^ Time instants can be inner oriented as either sources (see 
bottom left of the Inner orientation block in Figure]^ or sinks (see bootom right of the Inner orientation 
block in Figure]^. The outer orientation is induced by inner orientation, see the Outer orientation block in 
Figure For a more detailed discussion of orientation see [niEaiMi and especially the more recent work 
by Tonti m where the particular case of time is addressed. 
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Inner orientation 


Outer orientation 



Figure 2: Example of possible orientations of points and edges. See 11111521 l69l for a more detailed discussion of orientation. 


If we endow each of the geometric objects with a default orientation then we call the grid an oriented 
grid, see Figure for examples of how the geometric objects can be oriented. 

Given a grid, e.g. D, the space of /c-chains, Ck{D), is the collection of weighted fc-cells. A fc-chain, 
C(fe) G Ck{D), is a formal sum of fc-cells, T(^k),i G D, 

i 

By formal sum we mean a collection of cells and weights, {T(fcyi,c®}, where the d denote the weights. 

The boundary operator, d, on fc-chains, with fc > 1, is a homomorphism, d : Ck{D) —>■ Ck-i{D), such 
that: 


^^{k) ^ ^ ^ ^ ^ ^'^(k),i d ■ 

i i 

For 0-chains, 9c(o) = 0. The boundary of a /c-cell, is a (A: — l)-chain formed by its oriented faces. The 
coefficients of this (fc — I)-chain are associated to each of the faces and are given by the orientations: 

9T(^k),i = T(k-l),j j 

3 

with: 

{ el = I, if the orientation of T[k-i),j equals the default orientation 

el = —1, if the orientation of is opposite to the default orientation 

el = 0, if is not a face of j, 

and ddcf^k) = fl¬ 
it is possible to establish an isomorphism between the A:-chains and the corresponding column vectors of 
its coefficients and we use C(fc) to represent the column vector of the coefficients of the fc-chain C(fc), see [52] 
for more details. 

Example 2. As an example of k-chains and the boundary operator, consider the grid represented in Figure^ 
where the default orientation of the geometric objects is depicted. 





'^(O).O ''■(0),! ''’(0).2 


Figure 3: Example of a 1-dimensional grid, D, with the representation of the default orientation of the geometric objects. 
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Then we can construct the following 1-chains: 


C(l) = T(i),i +r(i),2 

“( 1 ) = 

and the following 0-chains: 

C(0) = '’■(0),0 + ■’■(O),! + '’■(0),2 


C(l) = 

«(!) = 


“(o) - - 




C(0) 


J( 0 ) 


1 

1 

0 

-1 


1 

1 

1 

0 

-1 

0 






T ( 0),0 


The boundary operator, d, applied to the 1-cochain, C(^i), results in: 

9c(i) = 


9c(i) =9r(i)_i + i9r(i),2 

= ~'’‘( 0).0 + ''■( 0 ),! — T -( 0),1 + T { 0),2 
= ~'''( 0),0 + ''■( 0),2 

and to the 1-cochain a(i) results in: 

aa(i) =-ar(i)_2 

= -(-T(o),i+'r(o),2) <— 

= '^(0),l ~ '’'(0),2 


-1 0 

1 -1 

0 1 


i9a(i) = 


-1 0 

1 -1 

0 1 


0 

-1 


- T ( 1),2 




-T( 0),1 


= E 


=( 0 . 1 ) 


(0,1) 


0 

-1 


T{0),2 


-1 

0 

1 


0 

1 

-1 


Remark 1. The matrices E(o i) present in Example^ are an example of incidence matrices, which 

are matrix representations of the boundary operator. 


For a grid to be considered a cell complex we impose that if a fc-dimensional geometric object, j, is 
in the grid, its boundary, dT(^k).ij is also in the grid. We see that the primal grid, D in Figure]^ is a cell 
complex. The dual grid, D in Figure is only a cell complex if boundary nodes are added. For a detailed 
discussion of this topic the reader is directed to [ISl [S 2 ] . 

Dual to the space of fc-chains, Ck{D), is the space of k-cochains, C^{D), defined as the set of homomor- 
phisms, : Cfe —>■ K, 


,(fc) 


,C(fc)) := (c(fe)) . 


With the duality pairing between cochains and chains, it is possible to define the formal adjoint of the 
boundary operator, the coboundary operator, 5 : C^{D) —>■ 


((5c('=),C(fe+i)) := (c('=),ac(fc+i)). 


(3.1) 


This equation is the discrete analogue of the generalized Stokes’ equation, (2.20). 

Analogous to the exterior derivative, the coboundary oper ator is nilpotent: 55c^^'> = 0 for all S 
C^{D). This property follows directly from its definition, (3.1), and from the nilpotency of the boundary 
operator := 0 for all S Ck{D). 

Let Ck{D) be the space of fc-chains with basis {T(k),j}, then the basis of the dual space of fc-cochains, 
Ck{D), is given by such that {T(k),j) = (here 5® are the coefficients of the Kronecker delta). 

In this way, all fc-cochains can be represented as linear combinations of these basis elements. 


cW = ^rW’®c,. 
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Similar to fc-chains, there exists an isomorphism between fc-cochains an the corresponding row vector of its 
coefhcients and we use to represent the row vector of the coefficients of the fc-cochain , see [52] for 
more details. With this isomorphism that identifies A:-chains and fc-cochains with their coefficients we can 
establish a matrix representation for the coboundary operator, 5: 

rank(C*’(D)) rank(C*;+i(£>)) 

(c«,9c(fc+i)) = ^ ^ c, (3.2) 

rank(C^(r>)) rank(Cfc+i(D)) 

= H (3.3) 

Z=1 j=l 

= (<5c«,c(fe+i)). (3.4) 

Therefore, when the row vector contains the expansion coefficients, c^, for the cochain then the 
row vector /j+i) contains the expansion coefficients for Sc^^\ 


4. Bridging the continuous and the discrete 

In Section [^ we have introduced the key concepts of differential geometry, integration over manifolds, 
differential operators and established fundamental integral relations. In Section[^we have presented an anal¬ 
ogous mathematical structure but at the discrete level, using algebraic topology. We now have two distinct 
worlds, one continuous and one discrete, with paralleled mathematical structures but still no connection 
between them. What we propose to do in this section is to explain how to convert continuous forms into 
discrete cochains and vice versa. In this way, we establish a bridge between the continuous and the discrete. 


4.I. From continuous to discrete 

The reduction operator, TZ : (Af) —t C^{D), maps fc-differential forms onto fc-cochains by 


TZa^'^^ 


:= (7^a«,T(fc),,) := / = (««,r(fc),i), & Ck{D). 


(4.1) 


Therefore, for all fc-chains, C(fc) € Ck{D), the reduction of a fc-form, € A^{Af), to a fc-cochain, 
aW g C'=(Z?), is 


,(fc) 


(c(fe)) := (7^a(''^C(fe)) = ^c*(7^a(''^r(fc),i) = = f 


(4.2) 


By defining the reduction operator in this manner an important commuting property relating the exterior 
derivative and the coboundary operator is satisfied: 


7^d = 5n. 


(4.3) 


Commuting relations are one of the key points of the mimetic framework, see [4bl I52j for a more ex¬ 
tensive discussion. In this case, this commuting relation states that taking the exterior derivative and then 
discretizing lead to the same result as discretizing first and then taking the discrete derivative. 

4-2. From discrete to continuous 

The operator that establishes the relation in the opposite direction, mapping fc-cochains into fc-differential 
forms, is the reconstruction operator, F : C^{D) —>■ where A^(A/') is the range of I. This operator 

can be implemented in different ways, see |S2] for a more detailed discussion and section Section 4.4 for its 
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use in this work. Despite the different possibilities to define this operator, the following commuting property 
must be satisfied, 

dZ = IS. (4.4) 


The reconstruction operator is essentially an interpolation operator which interpolates the discrete 
fc-cochains to continuous fc-forms. The commuting property (4.4) states that reconstructing and taking the 
exterior derivative is identical to first taking the discrete derivative and then reconstructing. This condition 
poses a severe restriction on the admissible reconstruction forms. For triangular elements these reconstruc¬ 
tion forms are known as Whitney forms, [13 [72]. The polynomial reconstruction forms on quadrilaterals 
are known as edge forms, 133112], and are the ones used in this work, see Section [4.4| 

Furthermore, the reconstruction operator, I, must be the right inverse of the reduction operator, 7^, so 
TZI = Id on C^{D) and it should approximate the left inverse of TZ, so ITZ = Id -I- 0{hP), where p is denotes 
the order of the approximation, for smooth differential forms. 

The first property, TZI = Id, is a consistency condition. The second property, ITZ = Id -|- 0{hP), is an 
approximability condition. 


4 . 3 . Mimetic projection 

The mimetic projection operator is obtained by combining the reduction and reconstruction operators: 


TTh :=IoTZ, 


(4.5) 


K\M) 

K 

C^iD) 



where A^(A/’; Ck) denotes the range of the reconstruction operator, I. 

Due to the commuting properties of the reduction and reconstruction operators we have: 

TThd = ITZd = ISTZ = dlTZ = di^h. 


Having defined the reduction, reconstruction andprojection operators for the primal grid D, equivalent 
operators can be defined for the dual grid, D, Figureld The space A^(A/’; Cu) is the space of discrete fc-forms, 
with fc-cochains associated with the dual fc-cells, 




Where Ck is the space of fc-chains on the dual complex associated to the dual grid, D, TZ is the reduction of 
differential forms on the dual chains and I constitutes the reconstruction of differential forms from cochains 
defined on the dual complex. 

There are essentially two different ways in which the discrete Hodge-* operator can be implemented. 
Either we can use the definition of reconstruction of differential forms or we can use the definition of the 
inner product of differential forms. 


For the first approach, we apply directly the definition of Hodge-* operator, (2.14) or its one dimensional 


version (2.16), to the reconstruction of a fc-cochain and subsequently reduce the result on the topological 
dual grid, i.e. TZ-kl or TZ-kl. Since here we use directly the reduction and reconstruction operators we refer 
to it as canonical Hodge. 


For the second approach we use (2.15): 


(»f,tf>) 


L2(aA) 


:= / a 
Jm 


(k) 


A */3 


(fc) 


and note that in this definition the Hodge-* is applied to the second argument, therefore the inner product 
induces the Hodge-* operator. Since the underlying principle of this construction of the Hodge-* operator 
is an inner product we refer to it as Galerkin Hodge. 

See Figure]^ for an example of the action of the two different Hodge-* operators on a 1-form. 
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Figure 4: Example of the action of the two different discrete Hodge-* operators. The two operators result in different ap¬ 
proximations (blue line and red line). If we refine the mesh or increase the polynomial order, the two operators will converge. 
Left: the analytical function and corresponding cochain (lines) associated to the degrees of freedom (dofs) of S A^CT ; Ci). 
Right: the reconstruction (lines) and corresponding cochains (dots) associated to the degrees of freedom of S A®(T;C'o)- 


4.4- Basis forms 

Consider a one dimensional domain, fl, and an associated grid consisting of a collection of points T(o)_i 
and line segments connecting the points T(i) i. 

Let be the space of smooth differentiable forms in il. Additionally, let the finite dimensional space of 
differential forms be defined as A^ = span }) ,1 = 1 ,..., dim (A^), where S A^ are basis fc-forms. 

Under these conditions it is possible, see (52], to define a projection operator tt/i which projects elements of 
A^ onto elements of A^ which satisfies, 

TThd = diTh. (4.6) 


Furthermore, it is possible to write. 


where 


TTha^^'> = ^ aie. 


(k) 

l^i 5 


f and f = i5*, A: = 0,1. 


(4.7) 


(4.8) 


4 . 4 . 1 . Primal grid 

Consider then a 0-form, (Qref), where Qref ■= f & [“I?!]- Define on Qref a cell com¬ 

plex D of order p consisting of [p + 1) nodes T(o),i = fi, where —l<fo<-"<ii<fp^l are the 
Gauss-Lobatto-Legendre (GLL) quadrature nodes, and p edges, T(i) j = ^i], of which the nodes are the 

boundaries. The projection operator, tt^ reads, 

^ 0^6-°^ (4.9) 

i=0 


where = k {f) are thep*^ degree Lagrange polynomials and ai = (fi). Note that (fj) = k {fj) = 

(4.8). In Figure|^ top left, an example of basis 0-forms of polynomial order four is presented. 
Similarly, for the projection of 1-forms in one dimension, [32 [62] derived 1-form edge polynomials, 

ef^ G AiiQref), 


Sj, as in 


ef^=e*(OdC, with e,(C) 


i-1 


E 


dlk 


i = l,...,p. 


(4.10) 
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note that we again have, 


= 






e, ( 0 d^ = <5;. 


(4.11) 


In Figure]^ top right, an example of the 1-form basis polynomial of order four is presented. 


4-4-^- Dual grid 

For the dual grid we again consider a 0-form defined on the reference interval, € hP{Qref)- In this 
case, on the reference interval, Qref, we define a cell complex D consisting of (p + 2) nodes ^(o),! = Cii where 

— 1 < Co < • • • < < Cp-i-i ^ {p+ 1) edges, = Cj-IjC* ; of which the nodes are the boundaries. 

The nodes used are the p nodes associated to Gauss-Legendre (GL) quadrature together with the two end 
points { — 1,1}. The projection operator, ixh reads. 


2 = 0 


(4.12) 


where = li (C) are the Lagrange polynomials of degree {p+ 1) and Ui = (^*)‘ ^ote that 

h (c,) ~ 


= (5!-, as in (4.8). In 


Figure 


bottom right, an example of basis 0-forms of polynomial order four is 


presented. 

Similarly, the dual I-form basis polynomials, € A^{Qref), are obtained from 
ef^=e*(C)dC, with = i = l,...,p+l, 


fc =0 


note that we again have, 


(4.13) 


(1) ^ 






e, (C)dC = <5;. 


In Figure]^ bottom left, an example of the I-form basis of polynomial order four is presented. 
4.4-3. Exterior derivative 

The exterior derivative of the primal basis 0-forms is given by. 


j (0) dli dlk dlk ( 


dC ^ dC 

fc =0 ^ fc =0 ^ 


( 1 ) , g(i) 

-1-1 w tj , 


f = 0, ...,p. 


(4.14) 


(4.15) 


The derivation of the exterior derivative of the dual basis 0-forms is obtained following the same procedure 
by substituting the primal functions by the respective dual functions. 

Example 3. Consider a 0-form, and the 1-form, = da^°\ both defined in = [—1,1] and take a 
grid D as the one defined in Figure If the discrete differential forms are given by: 






•- 

T (0).0 


''■( 0 ),! ''’( 0),2 

Figure 5: Representation of the 1-dimensional grid, Z), with the default orientation of the geometric objects. 
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then 


TTha^°^ = aoeQ°^ + + a 24 °^ + /3ie4\ 


d'KhOi^'^"’ = d ^aoEg^^ + aie[^^ + a 2 e 2 ^^J 
= agdeg^^ + aide^°^ + a2de2°^ 


= ao (-ei^^) + «i - 4^^) + «2 (e4^) 

= (ai - ao) + (02 - ai) 

= /3iel')+/3ie4) =/3(i). 

This means then that j3i = (ai — ao) and /32 = (a?. — ai). This corresponds exactly to the coefficients 
obtained from the projection of since, recall (j.l): 


■■= [ /?^^^ = / da(°) = a(T(o),*) - a(r(o),*_i), 


by Stokes’ theorem, (2.19). 

The exterior derivative of a discrete 0-form can then be written as, 


p p 


i—l j—Q 


( 1 ) 


(4.16) 


where, E^q are the coefficients of the matrix representation of the coboundary operator, (3.2), see 
for more details. 


5. Time integrators 
5.1. Introductory concepts 

As said before, the field of geometric integration is vast. Therefore, in what follows next, only key 
concepts essential to the understanding of the basic properties (such as symplecticity) of the numerical time 
integrators presented will be introduced. For a more detailed discussion and elaboration on these concepts, 
the authors suggest [371117]. 

Definition 1 (Symplectic form). |I]/ Let A4 be a manifold of dimension 2n. A symplectic form is a 
2 -form on At such that: 

1 . is closed, that is dw^^^ = 0; 

2. For each p S At and v,u £ TpM the identity uj^^\v,u) = 0 is satisfied if and only if v = 0 or u = 0 
(w^^^ is non-degenerate). 

For a Hamiltonian system, which involves a pair of n-dimensional variables q (configuration variables) 
and p (momenta), the symplectic form takes the form: 

= X! ^ 

i 

Definition 2 (Symplectic map). |IJ Given a 2n-dimensional manifold At, a diffeomorphism $ : 
At I —is a symplectic map if it preserves the symplectic form, that is: 

4>*a;(2) = 
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Figure 6: Basis functions associated to a p = 4 spectral finite element. Top left: primal grid 0-form basis functions (^). 
Top right: primal grid 1-form basis functions (^). Bottom left: dual grid 1-form basis functions (^). Bottom right: 
dual grid 0-form basis functions (^). 


Definition 3 (Hamiltonian vector field), On a 2n-dimensional manifold A4 a vector field vh ^ 
is called Hamoltonian if there is a 0-form € A^(At) such that: 

^ d^(0) ^ 

and is the symplectic form. Moreover, in local coordinates, vh takes the form: 


Vh 



dH d 
dpi dqi 


dH d \ 
dqi dpi J 


Additionally, it is possible to show, mm. that the integral curves generated by vh satisfy Hamilton’s 
equations: 

dqi dH dpi dH 

dt dpi ’ dt dqi ’ 


which is simply a restatement of (1.2). 
From these definitions it follows that: 


Theorem 1 (Properties of Hamiltonian systems). \37l Given a 2n-dimensional manifold j\4 and 
an Hamiltonian vector field vh associated to the Hamiltonian H^^\ then the flow generated by the vector 
field, has the following properties: 
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1 . is symplectic. 

2. The Hamiltonian is conserved along integral lines: . 

Proof 1. See |I]/. 


Hence, integral lines are isocontours of the Hamiltonian. Therefore, given the symplecticity and energy 
conservation properties of Hamiltonian systems it is natural to search for numerical integrators that are 
symplectic or exact energy preserving, leading to numerical integral lines that remain in a bounded region 
close to the real integral curve (symplectic integrator) or lie on it (energy conserving integrator). 

Definition 4 (Symplectic numerical one-step method). Jg7| / A one-step numerical method, <i);j : {qtg,pto) S 
I— {qt^,pti) G is symplectic if the one step map: 



is symplectic whenever the method is applied to a smooth Hamiltonian system. 


In this work we restrict our study to either symplectic or exact energy preserving time integrators. The 
choice between an integrator that is exact energy preserving or symplectic depends on the problem at hand, 
see for example [65j . 

The extension of the mimetic framework to time dependent problems, and in particular to the solution 
of first order ordinary differential equations, is done by noticing that the system of first order ordinary 


differential equations, (1.1): 


dy 

dt 


= Hy), y{to) = y°, yeMcW^ and tG[to,tf]c 


can be rewritten in differential geometric terms as: 


*dyW = y(°)(to) = 2/°, 2/f e AO([to,t/]), hf^ G A°([to,t/]) ■ 


(5.1) 


That is, each coordinate of an integral curve is a 0-form in time and its time derivative is a 0-form on the 
dual grid, obtained by exterior differentiation in time. 


Based on the formulation (5.1), two arbitrary order discrete integrators related to Gauss collocation 


methods can be constructed: one based on the canonical Hodge-* and another one based on the Galerkin 
Hodge-*. Each discrete integrator has different properties, as will be seen. 

Although we focus on Hamiltonian systems, the time integrators we will present can be applied to systems 
of any dimension (Hamiltonian and non-Hamiltonian) but, as will be shown, when applied to autonomous 
Hamiltonian systems, they display important properties. Therefore, for the sake of clarity, and without loss 
of generality, we present initially a first order ordinary differential equation only in one variable and then 
extend it to systems. For a system with one variable, the formulation presented in (|5.1[), becomes 


*dy(o) = hW(y(o)), G AO([to,t/]), h^°^ G A°{[to,tf]) . 


(5.2) 


5.2. Mimetic canonical integrator 

As was done in Section |4.4[ the physical quantities are discretized in the following way: 

y(o) g A0(r) ^ G A0(r) and G A°(r) ^ hf G Al{T) , 
where T = [tQ,tf]. The discrete physical quantities are represented as: 

3=0 3 = 1 


19 








where t = and e^°^(r), e^^\T) are the 0 -form basis functions on the primal and dual grids, 

, ..-(0) 

" we 


respectively, as described in Section 4.4 The only difference here being that for the dual grid basis. 


consider only the Gauss-Legendre nodes: 


and 


,-( 0 ) 


(T)=f®(T) j = l,...,p. 


Note also that: 


yj — y{o) and 






(5.4) 


ki=i 


where r-l are the (p -I- 1) intermediate time instants associated to the nodes of the primal grid, are the 
p intermediate time instants associated to the nodes of the dual grid and P = + 1) *^ 2 *° + ^o- It is 

important to remark that the upper index corresponds to degrees of freedom in time. A lower index will 
later correspond to distinct physical variables that appear in a system of ordinary differential equations. 
The choice to display an upper index instead of another lower index was done simply to make a distinction 
between space and time degrees of freedom. With this discretization the degrees of freedom appear staggered 
in time, as can be seen in Figure 



Figure 7: Geometric interpretation of the solution of |5.2i obtained with the mimetic discretization (5.101: (f, The 

nodes of the primal grid where the trajectory is discretized are represented in red (corresponding to Gauss-Lobatto-Legendre 
quadrature nodes). The nodes associated to the dual grid where the time derivative is discretized are depicted in blue (cor¬ 
responding to Gauss-Legendre quadrature nodes). The flow field, represented by arrows, is tangent to the curve at the Gauss 
nodes. That is, the derivative of the approximate trajectory is exactly equal to the flow field at the Gauss nodes. 


Taking the discretization presented in (5.3) and substituting in (5.2) we obtain: 

★ d 



(5.5) 


Given this discretization and using (4.161, it is possible to rewrite (5.5) as: 

p 


p 

^ E 

j = l,fc=0 


i=i 


(5.6) 


Following the canonical projection procedure as presented in Section 4.3 the Hodge-* of the term on the 
left can be computed, yielding: 


E 


j=l,/=l,/c=0 


zk,l 

-( 0 , 1)2 


ez(r^) ef\T) = ^ h^ef\T) 


i=i 


(5.7) 
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Once more, since the basis 0-forms are the same, this equality can be stated as: 


1^1,k^o 


(5.8) 


Explicitly expanding , (5.4), one obtains: 




(5.9) 


i=l,fc=0 


\.k=0 


This is a non-linear set of equations for the unknown values which can be solved, for instance, with a 
Newton-Raphson iterative solver, for known initial values = yitg). 

For a system of ordinary differential equations, as in (5.1), instead of a single time dependent unknown 
y{t), we now have M distinct time dependent unknowns yi, with i = 1,..., M. Following the same steps as 
presented previously, but now for multiple j/’s, one arrives at the following algebraic system: 

E E^o|i)?/fe;(f^) =f^ 2 /fe^°^(-f^),...,^?/^e[,°^(f^)V j = l,...,p and i = 

/^1,/c^O \k^0 k^O / 

(5.10) 

The numerical time integrator resulting from this discretization will be referred to as mimetic canonical 
integrator, or MCI, since it is obtained from a canonical projection 'Kh. 

Remark 2. It is important to note that this formulation, although different from Gauss-Legendre collocation 
methods, yields identical results (up to machine accuracy). Gauss-Legendre collocation method constructs an 
interpolating polynomial with Gauss-Legendre nodes for both the trajectory and its derivative, whereas here 
we use an interpolating polynomial with Gauss-Lobatto-Legendre nodes for the trajectory (primal grid) and 
another interpolating polynomial with Gauss-Legendre nodes for its derivative (dual grid). 

5.3. Mimetic Galerkin integrator 

As was done for the mimetic canonical integrator, the physical quantities are defined in the following 
way: 

y(0) g ^ yW ^ ^(0) g ^ 

where T = [tQ,tf]. The discrete physical quantities are represented as: 


yr = T.y^^T{r), 


(5.11) 


1=0 


where r = ^^= 3 ^ — 1 ^ and ej'^E) Ih® 0 -form basis functions on the primal grid, as described in 
Section 14.41 

= j= 0 ,...,P, . 


Note also that: 


= y^^1(t^) , 


(5.12) 


where t^ = (r^ -I- 1) *^ 2 *° + ^0 and are again the (p -I- 1) intermediate time instants associated with the 
nodes of the primal grid. 


For the Galerkin Hodge formulatio n, i nstead of using a canonical projection to g o fro m (5.6) to (5.7) 
and to discretize h^^\ one starts with (5.2), discretizes the left hand side and applies (2.15): 

E / efHT)A*e®(r) = > m = 1 ,... ,p , (5.13) 

1 — 1 1 -n P \ \ 1 -n / / ^ 


j = l,i=l,/c=0 
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where T = [—1,1]. Using (2.171 together with (2.16) the previous expression reduces to: 


H / Z|(T)r^(T)V5dT= / h\^y’^e^°\TUl<^^{T)y/^AT m = 1,... ,p , (5.14) 

j=i,;=i,fc=0 \fc=0 / 

where, for constant time steps, ^Jg = This is again a non-linear set of equations for the unknown 

values . 

As before, for a system of ordinary differential equations, as in (5.1), there exists M distinct time 
dependent unknowns yi, with i = 1,... ,M. Following the same steps as presented previously, but now for 
multiple y’s, one arrives at the following algebraic system: 

1 —1 u—n t/—1 JI \ —pj k—0 ' 




.... 

with TO = 1,... ,p and i = 1,..., M. The numerical time integrator resulting from this discretization will 
be referred to as mimetic Galerkin integrator, or MGI, since it is obtained from a Galerkin projection. 


Remark 3. It is important to note that this formulation turns out to be a reformulation of the energy 
preserving collocation method introduced in which in turn is an extension to higher order of the method 
introduced in m- 


5.4- Properties of the mimetic integrators 

Having presented the two mimetic integrators it is essential now to investigate their properties. In this 
section we will prove that the mimetic canonical integrator is symplectic and that the mimetic Galerkin 
integrator is energy preserving. 


5 . 4 . 1 . Symplecticity of mimetic canonical integrator 

We start by noting that the mimetic canonical integrator is a special case of a Runge-Kutta method. 


and let bi, aij (i, j = 1,..., s) be real numbers. An s-stage Runge-Kutta method is given by: 


Definition 5 (Runge-Kutta integrator), Take the system of ordinary differential 


equations in 



= y° -I At'^bik, . 
2=1 


i = l,...,s 


(5.16) 

(5.17) 


By setting ki = ^ (U), expressing ^ (t) and, by integration, y{t) in terms of ki it follows that the 
mimetic canonical integrator is a special case of a Runge-Kutta method with a full matrix of coefficients 
Qij, therefore an implicit Runge-Kutta method. We show this explicitly for a first order differential equation 
only in one variable. For that we start with (5.16) and replace the right hand side using (5.9): 


1^1,k^o 


(5.18) 


We can write this expression in matrix notation as: 


(5.19) 
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where = e\^\c)- 


This can be rewritten as: 


fe = A (Ey + cy°"j , 


(5.20) 


where, 


i = *=1. 

Ci — ^(0,1)’ * — 1) • ■ ■ jP ! 


f = y\ i = i,...,p. 


We can then rewrite (5.20) into: 


E'^A-ifc - E-icy° = y. 
Given the definition of E and c we have that: 


-E-1c2/0 = 


1 if i < j 
0 if i > j 
,0 


y 


From (5.24) it is straightforward to see that: 


yP = 


^(e^IA^i) fc,+ 2 / 0 . 


From (5.27) and (5.17) we can directly see that: 

bi = — 

At 


From (5.24) and (5.16) we have, 


“ At 


(e-^a-i) .. 

\ / p,l 

'' .7 


(5.21) 

(5.22) 

(5.23) 

(5.24) 

(5.25) 

(5.26) 


(5.27) 


Theorem 2 (Symplecticity of Runge-Kutta integrators), If a Runge-Kutta method conserves 
quadratic first integrals (i.e. /(j/i) = 7(2/o) whenever I{y) = y^Cy, with symmetric matrix C, is a first 
integral of (1.1^ ), then it is symplectic. 

Proof 2. See 


Remark 4. It is important to note that a quadratic first integral I{y) = y^Cy, with C constant, is conserved 

if: 

% ^ ti ^ ^ 2y*Ch(y) = 0, Vy . (5.28) 

We stress that this equality holds for all points in the phase space. 


Theorem 3 (Mimetic canonical integrator, (5.10), is symplectic). 

tor, (5.10), is symplectic. 


The mimetic canonical integra- 
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Proof 3. Consider a system of ordinary differential equations such as in (1.1) and take as the approxi¬ 
mated polynomial solution obtained by the mimetic canonical projection integrator. Assume I{y) = y*Cy is 
a quadratic invariant of the system (1.1), see \5.28 ). Therefore: 

Cf at r\n, 

I{y^{tf))-I{yM) = I = / 

The integral term on the right is a polynomial of degree at most (2p — 1), which can be exactly integrated 
by Gauss quadrature, over p collocation points (the Gauss nodes of the dual grid), t^, using the integration 
weights Wk. Therefore, one can write: 

I{y,{tj)) - I{y,ito)) = J2y^ (P)‘ (P) w, . 


By construction of the mimetic canonical integrator we have that: 

and consequently: 


I{yPtf)) - I{yPt,)) = ^ 2y, (P)‘ {P) wu C h {y, (P)) w, . 


Equation (5.28) states that if the continuous system of equations conserves a quadratic first integral then: 

yPh{y)=0, Vy. 

Therefore all the terms in the summation term on the right hand side are zero: 

Vh YY) =0- 


This directly implies that: 


Pyptf)) - liVkito)) = 0 ■ 


Therefore quadratic first integrals are conserved by the mimetic canonical projection integrator. Since it is 
a Runge-Kutta method, by Theorem\^it is symplectic. 

5 . 4 . 2 . Energy conservation of mimetic Galerkin integrator 

Consider the case of autonomous Hamiltonian systems, where y = {p,q), q & K"*, p G K™ and, 


dy 

dt 


= J VH{y), y{to) = y, yGMc 


D 2 m 


and t S / C M with J = 


0 U 

-Im 0 


(5.29) 


where \m is the mxm identity matrix. The solution obtained with the mimetic Galerkin integrator, (5.15) 


exactly preserves the energy of Hamiltonian systems. To prove that, we will follow a similar approach to 

[MIEt]. 


Theorem 4. The mimetic Galerkin integrator of (5.15) exactly preserves the energy of Hamiltonian sys¬ 
tems. 

Proof 4. From the fundamental theorem of calculus: 

H{yh{tf)) - H{yPto)) = [ dH{yPt)) , 
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where 'y{t) is the path along which the line integral is computed. Using the chain rule and expanding the 
terms on the right hand side this can written as: 


2m 

HiVhitf))-Hiy^ito)) =Y^ j 


fAv^,h\'dH 

Tr-(y/iW)dt, 


dt J dpi 


In differential form notation this can he written as: 


H{yh{tf)) - H[yh{h)) = A* ^ —(y^)J 


(5.30) 


where T = [—1,1]. Since: 


Yi ef^r) = i = l,...,2m, 

j=l,l=l,k=0 


(5.31) 


using (5.15) and noting that for Hamiltonian systems we have that = (j ^ and therefore, 

(0) 


2m 


dH 




U =1 


it is possible to write: 


E 


'-Y^y^eiir^) A ^ (J-^ VH)f A 


(5.32) 


where r = 1,... ,p and i = 1,..., 2m. We then notice that: 


N*., := / A*ef) = / ff (r)[f (r)V^dr . 


IT 


(5.33) 


This integral can be exactly evaluated by Gaussian quadrature. Since Ifir) is a polynomial of degree {p— 1), 
the product will be a polynomial of degree at most 2{p — 1). Using Gaussian quadrature over the p Gauss 
nodes associated to the 0-forms fi, the quadrature is exact for polynomials up to degree {2p — 1) and 
therefore exactly evaluates this integral. In this way, the matrix N becomes: 

ri _ _ ^ ^ 

Njj = J lf{T)Uj{T),/gdT = ^^f(T'')r®(f'')wfc = = SijWj . (5.34) 


k=l k=l 

. 1-1 


Which means that N zs a diagonal matrix and therefore is also a diagonal matrix whose elements 
are given by (N^^)^ . = SijwJ^. If we substitute this into (5.32), we get: 


ez(f^) = kr wk i (J-' A «(0) , 

i=l,fc=0 


(5.35) 


r—1 


tr 


with j = 1 ,... ,p and i = 1,..., 2m. 


Now, substituting (5.35) in (5.31) yields: 

p / p c 


E ( E ^iT ^ f A *4°^ I = *kh,i: i = l,...,2m. 

j=l \r=l 


(5.36) 
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A*(yH)[°'> . (5.37) 


Substituting this expression into the right hand side of {5.30) yields: 

2m „ p / r \ 

H{y^{tf))-H{y^(t^))= E / E [SjrWf^ A*e^°A ef 

2_1 _1 ^ j _1?^_ ^ 

Rearranging the terms and summing in j, the following expression is obtained: 

H{y^{tf))-H{yM)= E E f / (Vi7)i°^ A *6®) J.fc f / A * (Vi7)f)) . (5.38) 

i=l,fe=lm=l / \Jt / 

The right hand side of this expression is equal to zero due to the skew-symmetry of matrix J, this equality 
becomes: 

HiVhitf)) - Hiy^ifo)) = 0 , 

and consequently the mimetic Galerkin integrator is exactly conserves the energy of the Hamiltonian system. 


Definition 6 (Symmetric integrator). Take the system of ordinary differential equations in (1.1). The 
numerical method <i>At : f/o G K" yi G M", which approximates this system of equations, maps points 
yQ in the instant to to points yi in the time instant ti = to At. The numerical method is symmetric if 
$_At o $At = I■ 


One can see that the canonical and Galerkin integrators are symmetric. For that recall (5.9) and (5.14) and 


note that substituting yh{—t) in (5.9) and (5.14) it solves both systems of equations for — At. 


6. Numerical tests 

In order to illustrate the properties of the mimetic integrators just presented we apply them to the 
solution of four test cases. We have chosen the circular motion problem (a Hamiltonian linear problem), 
the Lotka-Volterra problem (a non-Hamiltonian, non-linear polynomial problem), the pendulum problem (a 
Hamiltonian, non-linear and non-polynomial problem) and finally the two-body Kepler problem (a higher 
order Hamiltonian problem). 

6.1. Circular motion problem 

The mimetic integrators were applied to a circular motion problem, a linear problem. 


dP = o 
di y 

dt A 


( 6 . 1 ) 


whose Hamiltonian is given by H{p, q) = \ + P^) ■ 

In Figure 1^ one can see the numerical solution of the circular motion problem, (6.1), using the two 
mimetic integrators, with polynomial order in time pt = 2, a Runge-Kutta of 4th order and the simple 
explicit Euler. From this figure one can see that both mimetic integrators are able to perform a long time 
simulation, as they retain the behavior of the analytical solution. On the other hand, as expected, the explicit 
Euler scheme rapidly spirals outwards and the Runge-Kutta of 4th order quickly spirals inwards reaching 
the stable point (0,0). This is an important example that stresses the relevance of symplectic and energy 
preserving integrators. Although the Runge-Kutta integrator has the same convergence rate of the mimetic 
integrators, see Figure]^ (right), it completely fails to capture the circular motion of this system of equations. 
Both mimetic integrators exactly (up to machine precision) preserve the constant of motion E? = q^ -\-p^. 
In the case of the mimetic canonical integrator this happens because this quantity is a quadratic expression 
of the variables q and p, and this integrator is symplectic. For the mimetic Galerkin integrator this happens 
because this quantity differs from the Hamiltonian of the system by only a multiplicative constant. An 
interesting point to note is that, in this case, due to the linearity of the system, both mimetic integrators 
yield the same result, up to machine accuracy, as can be seen in Figure 10 Another relevant point to note 
is that, although the trajectory computed from both mimetic integrators is of order pt = 2 they show a 
convergence rate of 2p(. 
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Figure 8: Phase-space numerical (in color) and analytical (in black) solutions of circular motion, l |6.1| l, with initial condition 
yi = 2, y 2 = 0 and time step At = Is. Top left: mimetic canonical integrator with poynomial order in time pt = 2 (MCI2). 
Top right: mimetic Galerkin integrator with polynomial order in time pt = 2 (MGI2). Genter left: explicit Euler scheme, which 
clearly diverges. Center right: symplectic Euler scheme, which clearly diverges. Bottom: explicit Runge-Kutta of 4th order, 
which also clearly diverges. 
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Figure 9: Left: Error in time of R = for the numerical solution of the circular motion, (6.11, for At = Is, with the 

mimetic canonical integrator of polynomial order in time pt = 2 (MCI2), mimetic Galerkin integrator of polynomial order in 
time Pt = 2 (MGI2), Runge-Kutta of 4th order, explicit Euler scheme and symplectic Euler method. Since the quantity R is 
a function (square root) of a quadratic term, the mimetic canonical integrator exactly (up to machine accuracy) preserves this 
quantity along the trajectory, as can be seen. The same can be seen for the mimetic Galerkin integrator, since ^ {y^ + y^) is 
also the Hamiltonian of the system. Right: Gonvergence of the error of the path in phase space as a function of the time step 
(^) for mimetic canonical integrator of polynomial order in time pt = 2 (MCI2), mimetic Galerkin integrator of polynomial 
order in time pt = 2 (MGI2), Runge-Kutta of 4th order, explicit Euler and symplectic Euler. It is possible to observe the 4th 
order rate of convergence of the mimetic integrators. 
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Figure 10: Difference between the numerical solution of the circular motion, ( |6.1| l, for At = Is, obtained with the mimetic 
canonical integrator of polynomial order in time pt = 2 (MGI2) and the mimetic Galerkin integrator of polynomial order in 
time Pt = 2 (MGI2). For linear systems of ordinary differential equations, both solutions are equal up to machine accuracy. 
In this particular case, both integrators give the same result. This is due to the fact that the in system of equations ( |6.1[ | is 
linear. 
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6.2. Lotka-Volterra problem 

The Lotka-Volterra problem is governed by the following equations, 

^ = yi{y2-2) 

^ = 2/2(1-yi) ’ 


( 6 . 2 ) 


a non-linear polynomial problem, which is not a Hamiltonian system but conserves the quantity V(jji,y 2 ) = 
— 2/1 -|- log 2/1 — 2/2 + 2 log 2/2 along its trajectories. 


The Lotka-Volterra system of equations, (6.21, is an interesting problem to analyze since it is neither 


Hamiltonian (it is a Poisson system) nor has a conserved quadratic quantity. Here one can see that the 
trajectories obtained with both mimetic integrators and the Runge-Kutta of 4th order are very similar, see 
and the error as a function of time of the conserved quantity V = — 2/1 + log 2/1 —2/2+2 log 2/2, see 
(left). A similar behaviour is seen for the convergence plots as a function of time step size, see 


11 


12 


Figure 
Figure 

Figure 12 (right). In this case the advantage of using the mimetic integrators when compared to Runge-Kutta 
of 4th order does not show itself. 


29 







a Euler 
-Analytical 



2 

Vi 




Figure 11: Phase-space numerical (in color) and analytical (in black) solutions of Lotka-Volterra system, l |6.2| | with initial 
condition = 3, y 2 = 3 and time step At = 0.3s. Top left: mimetic canonical integrator with polynomial order in time pt = 2 
(MCI2). Top right: mimetic Galerkin integrator with polynomial order in time pt = 2 (MGI2). Center left: explicit Euler 
scheme, where one can clearly see the fast divergence of the solution. Center right: symplectic Euler scheme, where one can 
also see the divergence of the solution. Bottom: explicit Runge-Kutta of 4th order. 
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Figure 12: Left: Error in time of the conserved quantity V = —y\ + logt/i — 2/2 + 21ogy2 for the numerical solution of the 
Lotka-Volterra system of equations, ( |6.2[ |, for At = 0.3s, with the mimetic canonical integrator of polynomial order in time 
pt = 2 (MCI2), mimetic Galerkin integrator of polynomial order in time pt = 2 (MGI2), Runge-Kutta of 4th order, explicit 
Euler scheme and symplectic Euler method. Since this quantity V is not a quadratic term, the mimetic canonical integrator is 
unable to exactly preserve this quantity along the trajectory, as can be seen by the increase in the error with time. The same 
can be seen for the mimetic Galerkin projection integrator, since V is not the Hamiltonian of the system. Right: Convergence 
of the error of the path in phase space as a function of ^ for mimetic canonical integrator of polynomial order in time pt = 2 
(MCI2), mimetic Galerkin integrator of polynomial order in time pt = 2 (MGI2), Runge-Kutta of 4th order, explicit Euler and 
symplectic Euler. It is possible to observe the 4th order rate of convergence of the mimetic integrators. 
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6.3. Pendulum problem 

The pendulum problem is a non-linear, non-polynomial problem, 

/ ^ = -losing 

I da _ „ 


(6.3) 


with Hamiltonian H{p,q) = | — 1 0 cos g. 

The pendulum problem, (6.31 is a Hamiltonian system and here both mimetic integrators show again 
a clear advantage when compared to the Runge-Kutta integrator. As can be seen in Figure |13[ both 
mimetic integrators are capable of keeping the numerical trajectories close to the analytical solution, therefore 
retaining the qualitative behavior of the problem, whilst the Runge-Kutta integrator shows an inward 
spiraling behavior. This numerical solution is characteristic of a damped pendulum problem (with friction) 
but the problem being solved is a frictionless pendulum. This issue is relevant since the numerical method 
interferes with the physics of the problem being solved. 


32 






Figure 13: Phase-space numerical (in color) and analytical (in black) solutions of pendulum system, | |6.3[ l, with initial condition 
q = p = 0 and time step At = 0.4s. Top left: mimetic canonical integrator with polynomial order in time pt = 2 (MCI2). 
Top right: mimetic Galerkin integrator with polynomial order in time pt = 2 (MGI2). Genter left: explicit Euler scheme, which 
clearly diverges, with the system rapidly gaining energy. Center right: symplectic Euler scheme, which keeps a stable orbit but 
with a large difference in comparison to the analytical trajectory. Bottom: explicit Runge-Kutta of 4th order, which clearly is 
dissipative. 
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Figure 14: Left: Error in time in the Hamiltonian H = —10 cos q for the numerical solution of the pendulum problem, (6.3 i, 

for At = 0.1s, with the mimetic canonical integrator of polynomial order in time pt = 2 (MCI2), mimetic Galerkin integrator 
of polynomial order in time pt = 2 (MGI2), Runge-Kutta of 4th order, explicit Euler scheme and symplectic Euler method. 
Since the Hamiltonian, H, is not a quadratic term, the mimetic canonical integrator cannot exactly preserve this quantity along 
the trajectory, as can be seen. Nevertheless, the error remains bounded. On the other hand the mimetic Galerkin integrator 
exactly (up to machine precision) preserves it, since it is an exact energy preserving integrator. Both the explicit Euler and 
the Runge-Kutta methods, show an increase of the error with time. The symplectic Euler method, shows also that the error 
in the energy of the system is kept bounded. Right: Gonvergence of the error in as a function of ^ for mimetic canonical 
integrator of polynomial order in time pt = 2 (MCI2), mimetic Galerkin integrator of polynomial order in time pt = 2 (MGI2), 
Runge-Kutta of 4th order, explicit Euler and symplectic Euler. It is possible to observe the 4th order rate of convergence of the 
mimetic canonical integrator and the exact solution for the mimetic Galerkin integrator. The Runge-Kutta integrator shows a 
5th order rate of convergence. 
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6.4- Two-body Kepler problem 

Finally, the mimetic time integrators were applied to the solution of the two-body Kepler problem, 


dt 

dt 

dpi _ _ 
dt 

dp2 _ _ 
dt 


qi 


Q2 


(6.4) 


a higher order Hamiltonian system with Hamiltonian H{pi,p 2 ,qi,q 2 ) = 5 {pi + P 2 ) - /I 2 • 

This problem is also a Hamiltonian system and again both mimetic integrators show a clear advantage 
when compared to the Runge-Kutta integrator. As can be seen in Figure [Tsl both mimetic integrators are 
capable of keeping the numerical trajectories close to the analytical solution, therefore retaining the quali¬ 
tative behavior of the problem, whilst the Runge-Kutta integrator shows a slow inward spiraling behavior, 
once more indicating an artificial damping effect. The symplectic Euler shows a precession in the trajectory. 


Looking at the error in the Hamiltonian with time,Figure 16 (left), one can clearly see that the mimetic 
Galerkin integrator exactly preserves, up to machine accuracy, the Hamiltonian. The mimetic canonical 
integrator, although not preserving the Hamiltonian exactly, keeps its error bounded. On the other hand, 
the Runge-Kutta and the Euler integrators quickly build up the error in the Hamiltonian. In Figure 16 
(right) one can see the convergence rates of the different integrators. The Runge-Kutta shows a convergence 
rate close to order five and the mimetic canonical integrator shows a convergence rate of order four. In both 
cases the mimetic Galerkin integrator exactly conserves the Hamiltonian, therefore the convergence rate is 
shown as zero since the error is always zero. 
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Figure 15: Phase-space numerical solution (in color) of Kepler two body system, ( |6.4[ | with initial condition qi = 0.4, pi = 0, 
^2 = 0 and p 2 = 2 and time step At = 0.1s. Top left: mimetic canonical integrator with polynomial order in time pt = 2 
(MCI2). Top right: mimetic Galerkin integrator with polynomial order in time pt = 2 (MGI2). Genter left: explicit Euler 
scheme, where clearly the trajectory diverges, with the system rapidly gaining energy. Center right: symplectic Euler scheme, 
where one can see the precession in the trajectory. Bottom: explicit Runge-Kutta of 4th order, which diverges, loosing energy. 
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Figure 16: Left: Error in time in the Hamiltonian H = ^ (p| +P 2 )- , ^ for the numerical solution of the Kepler two 

V 91+92 

body problem, | |6.4| l, for At = 0.1s with the mimetic canonical integrator of polynomial order in time pt = 2 (MCI2), mimetic 
Galerkin integrator of polynomial order in time pt = 2 (MGI2), Runge-Kutta of 4th order, explicit Euler scheme and symplectic 
Euler method. Since the Hamiltonian, is not a quadratic term, the mimetic canonical integrator cannot exactly preserve 
this quantity along the trajectory, as can be seen. Nevertheless, the error remains bounded. On the other hand the mimetic 
Galerkin integrator exactly (up to machine precision) preserves it, since it is an exact energy preserving integrator. Both the 
explicit Euler and the Runge-Kutta methods, show an increase of the error with time. The symplectic Euler manages to keep 
the error in the energy bounded. Right: Gonvergence of the error in as a function of the number of time steps for mimetic 
canonical integrator of polynomial order in time pt = 2 (MGI2), mimetic Galerkin integrator of polynomial order in time 
Pt = 2 (MGI2), Runge-Kutta of 4th order, explicit Euler and symplectic Euler. It is possible to observe the 4th order rate of 
convergence of the mimetic canonical integrator and the exact solution for the mimetic Galerkin integrator. The Runge-Kutta 
integrator shows a 5th order rate of convergence. 
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7. Summary and outlook 

In this work we took as a starting point the mimetic framework established in [461152] and applied it to 
the numerical solution of systems of ordinary differential equations. This framework is based on two well 
established mathematical constructs: differential geometry (for the continuous formulation) and algebraic 
topology (for the discrete formulation). The connection between the two is established by a projection 
operator. With this framework, the topological relations can be exactly represented at the discrete level 
and all approximations lie in the metric dependent operators (e.g. Hodge-*). Moreover, these metric 
dependent operators do not have a unique discrete representation. Therefore, for example, depending on 
the representation of the Hodge-* operator, it is possible to construct discretizations with different properties. 

When this framework is applied to the solution of systems of ordinary differential equations, we have 
shown that two different time integrators can be obtained, one symplectic and the other exact energy 
preserving. The two integrators are obtained simply by choosing a specific discrete representation of the 
Hodge-* operator. For a discrete Hodge-* based on a canonical Hodge the integrator is symplectic, for a 
discrete Hodge-* based on a Galerkin Hodge the integrator is exact energy preserving. 

In future work, we intend to build a physically accurate framework for time dependent systems of partial 
differential equations. Furthermore, the symmetry property that characterizes both integrators presented 
can play an important role, for instance, in the application of the mimetic framework to the advection 
equation. Moreover, given the properties of the integrators presented in this work, and the properties 
obtained for pure spatial discretizations, [451152] , we also intend to explore a space-time formulation of the 
mimetic framework. 
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